Theme Song


1 SCSS Setup

# install.packages("remotes")
require('BBmisc')
#remotes::install_github("rstudio/sass")
lib('sass')
## sass 
## TRUE
/* https://stackoverflow.com/a/66029010/3806250 */
h1 { color: #002C54; }
h2 { color: #2F496E; }
h3 { color: #375E97; }

/* ----------------------------------------------------------------- */
/* https://gist.github.com/himynameisdave/c7a7ed14500d29e58149#file-broken-gradient-animation-less */
.hover01 {
  /* color: #FFD64D; */
  background: linear-gradient(155deg, #EDAE01 0%, #FFEB94 100%);
  transition: all 0.45s;
  &:hover{
    background: linear-gradient(155deg, #EDAE01 20%, #FFEB94 80%);
    }
  }

.hover02 {
  color: #FFD64D;
  background: linear-gradient(155deg, #002C54 0%, #4CB5F5 100%);
  transition: all 0.45s;
  &:hover{
    background: linear-gradient(155deg, #002C54 20%, #4CB5F5 80%);
    }
  }

.hover03 {
  color: #FFD64D;
  background: linear-gradient(155deg, #A10115 0%, #FF3C5C 100%);
  transition: all 0.45s;
  &:hover{
    background: linear-gradient(155deg, #A10115 20%, #FF3C5C 80%);
    }
  }
## https://stackoverflow.com/a/36846793/3806250
options(width = 999)
knitr::opts_chunk$set(class.source = 'hover01', class.output = 'hover02', class.error = 'hover03')



2 受講生によるテスト：MCMC algorithms and density estimation

課題をすぐに提出してください

課題の提出期限は、5月31日 15:59 JSTですが、可能であれば1日か2日早く提出してください。 早い段階で提出すると、他の受講生のレビューを時間内に得る可能性が高くなります。



2.1 説明

The R dataset faithful contains data on waiting time between eruptions (the column named waiting) and the duration of the eruption (the column named eruptions) for the famous Old Faithful geyser in Yellowstone National Park, Wyoming, USA.

In this case, you are asked to modify the MCMC algorithm provided in “Sample code for density estimation problems” (as opposed to the EM algorithm you used in the previous peer assignment) to provide a (Bayesian) density estimate the marginal distribution of the duration of the eruptions using a location-and-scale mixture of 2 univariate Gaussian distributions (as opposed to the location mixture of 6 univariate Gaussian distributions that we used for the galaxies dataset). Assume that the priors are \(\omega∼Beta(1,1)\), \(\mu_k∼Normal(η,\tau^2)\) and \(1/σ^2_k∼Gamma(d,q)\), where \(\eta\), \(\tau^2\), \(d\) and \(q\) are selected using an empirical Bayes approach similar to the one we used in “Sample code for density estimation problems”.



2.1.1 Review criteria

Reviewers will check whether the code has been modified correctly, and whether the density estimate you generate appears correct. Please remember that you are being asked to use a location-and-scale mixture to generate the density estimate, so the “Sample code for density estimation problems” cannot be used directly and requires some modification. Before submitting your answer, it might be useful to compare the density estimate generated by your algorithm against a kernel density estimate generated by the R function density(), and agains the answer to the previous peer assignment.



2.2 自分の提出物

2.2.1 Setup

if(!suppressPackageStartupMessages(require('BBmisc'))) {
  install.packages('BBmisc', dependencies = TRUE, INSTALL_opts = '--no-lock')
}
suppressPackageStartupMessages(require('BBmisc'))
# suppressPackageStartupMessages(require('rmsfuns'))

pkgs <- c('devtools', 'knitr', 'kableExtra', 'tidyr', 
          'readr', 'lubridate', 'data.table', 'reprex', 
          'timetk', 'plyr', 'dplyr', 'stringr', 'magrittr', 
          'tdplyr', 'tidyverse', 'formattable', 
          'echarts4r', 'paletteer')

suppressAll(lib(pkgs))
# load_pkg(pkgs)

## Set the timezone but not change the datetime
Sys.setenv(TZ = 'Asia/Tokyo')
## options(knitr.table.format = 'html') will set all kableExtra tables to be 'html', otherwise need to set the parameter on every single table.
options(warn = -1, knitr.table.format = 'html')#, digits.secs = 6)

## https://stackoverflow.com/questions/39417003/long-vectors-not-supported-yet-abnor-in-rmd-but-not-in-r-script
knitr::opts_chunk$set(message = FALSE, warning = FALSE)#, 
                      #cache = TRUE, cache.lazy = FALSE)

rm(pkgs)
dat <- fread('data/fuses.csv') %>% as.matrix
head(dat)
##            V1
## [1,] 1.062163
## [2,] 1.284631
## [3,] 2.733352
## [4,] 2.284062
## [5,] 3.289444
## [6,] 3.113181

Source : 400 x 1

2.2.2 Assignment

2.2.3 Marking

Is the sampler for \(c\) correct?

The full conditional for \(c_i\) is given by,

\[\begin{align} \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \end{align}\]

(This assumes that the mixture weights are parameterized as

for(i in 1:n){
  v    = rep(0,2)
  v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
  v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
  v    = exp(v - max(v))/sum(exp(v - max(v)))
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
## Error in 1:n: NA/NaN argument
  • 0点 No
  • 1点 Yes

Is the sampler for \(\omega\) correct?

The full conditional is simply

\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]

The following one line of code implements it:

w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
## Error in rbeta(1, aa[1] + sum(cc == 1), aa[2] + sum(cc == 2)): object 'aa' not found
  • 0点 No
  • 1点 Yes

Are the sampler for \(\mu\) correct?

\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]

The following code implements the sampler:

# Sample the means
for(k in 1:2){
  nk       = sum(cc==k)
  xsumk    = sum(x[cc==k])
  tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
  mu.hat   = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
  mu[k]    = rnorm(1, mu.hat, sqrt(tau2.hat))
}
## Error in eval(expr, envir, enclos): object 'cc' not found

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Are the samplers for \(\sigma_k\) correct?

The full conditional takes the form:

\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]

The following code implements the sampler:

# Sample the variances
for(k in 1:2){
  nk    = sum(cc==k)
  dd.star = dd + nk/2
  qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
  sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}
## Error in eval(expr, envir, enclos): object 'cc' not found

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Does the code to generate the density estimate appear correct?

This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
## Error in apply(density.mcmc, 2, quantile, 0.025): object 'density.mcmc' not found
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
## Error in apply(density.mcmc, 2, quantile, 0.975): object 'density.mcmc' not found
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
## Error in plot(xx, density.mcmc.m, type = "n", ylim = c(0, max(density.mcmc.uq)), : object 'xx' not found
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
## Error in xy.coords(x, y, setLab = FALSE): object 'density.mcmc.lq' not found
lines(xx, density.mcmc.m, col="black", lwd=2)
## Error in lines(xx, density.mcmc.m, col = "black", lwd = 2): object 'xx' not found
points(x, rep(0,n))
## Error in points(x, rep(0, n)): object 'x' not found

Note that this will generate not only a point estimate, but also pointwise credible intervals.

  • 0点 No
  • 1点 Yes

Does the graph appear to be correct?

The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.

The full code for the MCMC algorithm follows:

x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))

## Priors set up using an "empirical Bayes" approach
aa  = rep(1,2)  
eta = mean(x)    
tau = sqrt(var(x))
dd  = 2
qq  = var(x)/2

## Initialize the parameters
w     = 1/2
mu    = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc    = sample(1:2, n, replace=T, prob=c(1/2,1/2))

## Number of iterations of the sampler
rrr   = 12000
burn  = 2000

## Storing the samples
cc.out    = array(0, dim=c(rrr, n))
w.out     = rep(0, rrr)
mu.out    = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost   = rep(0, rrr)

for(s in 1:rrr){
  # Sample the indicators
  for(i in 1:n){
    v    = rep(0,2)
    v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
    v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
    v    = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
  }
  
  # Sample the weights
  w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  
  # Sample the parameters of the components
  for(k in 1:2){
    # Sample the means
    nk    = sum(cc==k)
    xsumk = sum(x[cc==k])
    tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
    mu.hat  = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
    mu[k]   = rnorm(1, mu.hat, sqrt(tau2.hat))
    # Sample the variances
    dd.star = dd + nk/2
    qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
    sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
  }
  
  # Store samples
  cc.out[s,]    = cc
  w.out[s]     = w
  mu.out[s,]    = mu
  sigma.out[s,] = sigma
  logpost[s] = 0
  for(i in 1:n){
    if(cc[i] == 1){
      logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
    }else{
      logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
    }
  }
  logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
  for(k in 1:2){
    logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
  }
  if(s/500==floor(s/500)){
    print(paste("s =",s))
  }
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx  = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
  density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) + 
                               (1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

  • 0点 No
  • 1点 Yes



2.3 ピアレビュー

2.3.1 1st Peer

2.3.1.1 Assignment

prompt caption-text color-secondary-text

### Get a "Bayesian" kernel density estimator based on the same location mixture of 6 normals
## Priors set up using an "empirical Bayes" approach
aa      = rep(1,KK)  
## Error in eval(expr, envir, enclos): object 'KK' not found
eta     = mean(x)    
tau     = sqrt(var(x))
dd      = 2
qq      = var(x)/KK
## Error in eval(expr, envir, enclos): object 'KK' not found
mu_0    = rnorm(KK, eta, tau)
## Error in rnorm(KK, eta, tau): object 'KK' not found
sigma_0 = rinvgamma(KK, dd, qq)
## Error in rinvgamma(KK, dd, qq): could not find function "rinvgamma"
## Initialize the parameters
w     = rep(1,KK)/KK
## Error in eval(expr, envir, enclos): object 'KK' not found
mu    = rnorm(KK, mean(x), sd(x))
## Error in rnorm(KK, mean(x), sd(x)): object 'KK' not found
sigma = sd(x)/KK
## Error in eval(expr, envir, enclos): object 'KK' not found
cc    = sample(1:KK, n, replace=T, prob=w)
## Error in sample(1:KK, n, replace = T, prob = w): object 'KK' not found
## Number of iterations of the sampler
rrr   = 12000
burn  = 2000

## Storing the samples
cc.out    = array(0, dim=c(rrr, n))
w.out     = array(0, dim=c(rrr, KK))
## Error in array(0, dim = c(rrr, KK)): object 'KK' not found
mu.out    = array(0, dim=c(rrr, KK))
## Error in array(0, dim = c(rrr, KK)): object 'KK' not found
sigma.out = rep(0, rrr)
logpost   = rep(0, rrr)
for(s in 1:rrr){
  
  # Sample the indicators
  for(i in 1:n){
    v = rep(0, KK)
    
    for(k in 1:KK){
      v[k] = log(w[k]) + dnorm(x[i], mu[k], sigma, log = TRUE)  #Compute the log of the weights
    }
    v      = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i]  = sample(1:KK, 1, replace = TRUE, prob = v)
  }
  
  # Sample the weights
  w = as.vector(rdirichlet(1, aa + tabulate(cc, nbins=KK)))
  
  # Sample the means
  for(k in 1:KK){
    nk       = sum(cc==k)
    xsumk    = sum(x[cc==k])
    tau2.hat = 1/(nk/sigma^2 + 1/sigma_0[k]^2)
    mu.hat   = tau2.hat*(xsumk/sigma^2 + mu_0[k]/sigma_0[k]^2)
    mu[k]    = rnorm(1, mu.hat, sqrt(tau2.hat))
  }
  
  # Sample the variances
  dd.star = dd + n/2
  qq.star = qq + sum((x - mu[cc])^2)/2
  sigma   = sqrt(1/rgamma(1, dd.star, qq.star))
  
  # Store samples
  cc.out[s,]   = cc
  w.out[s,]    = w
  mu.out[s,]   = mu
  sigma.out[s] = sigma
  
  for(i in 1:n){
    logpost[s] = logpost[s] + log(w[cc[i]]) + dnorm(x[i], mu[cc[i]], sigma, log = TRUE)
  }
  logpost[s] = logpost[s] + log(ddirichlet(w, aa))
  for(k in 1:KK){
    logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = TRUE)
  }
  logpost[s] = logpost[s] + dgamma(1/sigma^2, dd, qq, log = TRUE) - 4*log(sigma)
  
  if(s/500==floor(s/500)){
    print(paste("s =",s))
  }
}
## Error in eval(expr, envir, enclos): object 'KK' not found

Provide code to generate the density estimate on a grid consisting of 150 points in the interval [0,7].

xx  = seq(0,7,length=150)
nxx = length(xx)

## Compute the samples of the density over a dense grid
density.mcmc = array(0, dim = c(rrr-burn, length(xx)))
for(s in 1:(rrr-burn)){
    for(k in 1:KK){
        density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn,k]*dnorm(xx,mu.out[s+burn,k],sigma.out[s+burn])
    }
}
## Error in eval(expr, envir, enclos): object 'KK' not found
density.mcmc.m = apply(density.mcmc , 2, mean)

Provide the a graph of the density estimate on the interval [0,7].

## Plot Bayesian estimate with pointwise credible bands along with kernel density estimate and frequentist point estimate
colscale = c("black", "blue", "red")
yy = density(x)
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)

plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruption Duration", ylab="Density")

polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col=colscale[1], lwd=2)
#lines(xx, density.EM, col=colscale[2], lty=2, lwd=2)
lines(yy, col=colscale[3], lty=3, lwd=2)
points(x, rep(0,n))
legend(5, 0.45, c("KDE","EM","MCMC"), col=colscale[c(3,2,1)], lty=c(3,2,1), lwd=2, bty="n")

2.3.1.2 Marking

Is the sampler for \(c\) correct?

The full conditional for \(c_i\) is given by,

\[ \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \]

(This assumes that the mixture weights are parameterized as

for(i in 1:n){
  v    = rep(0,2)
  v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
  v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
  v    = exp(v - max(v))/sum(exp(v - max(v)))
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
  • 0点 No
  • 1点 Yes

Is the sampler for \(\omega\) correct?

The full conditional is simply

\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]

The following one line of code implements it:

w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  • 0点 No
  • 1点 Yes

Are the sampler for \(\mu\) correct?

\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]

The following code implements the sampler:

# Sample the means
for(k in 1:2){
  nk       = sum(cc==k)
  xsumk    = sum(x[cc==k])
  tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
  mu.hat   = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
  mu[k]    = rnorm(1, mu.hat, sqrt(tau2.hat))
}

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Are the samplers for \(\sigma_k\) correct?

The full conditional takes the form:

\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]

The following code implements the sampler:

# Sample the variances
for(k in 1:2){
  nk    = sum(cc==k)
  dd.star = dd + nk/2
  qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
  sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Does the code to generate the density estimate appear correct?

This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

Note that this will generate not only a point estimate, but also pointwise credible intervals.

  • 0点 No
  • 1点 Yes

Does the graph appear to be correct?

The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.

The full code for the MCMC algorithm follows:

x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))

## Priors set up using an "empirical Bayes" approach
aa  = rep(1,2)  
eta = mean(x)    
tau = sqrt(var(x))
dd  = 2
qq  = var(x)/2

## Initialize the parameters
w     = 1/2
mu    = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc    = sample(1:2, n, replace=T, prob=c(1/2,1/2))

## Number of iterations of the sampler
rrr   = 12000
burn  = 2000

## Storing the samples
cc.out    = array(0, dim=c(rrr, n))
w.out     = rep(0, rrr)
mu.out    = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost   = rep(0, rrr)

for(s in 1:rrr){
  # Sample the indicators
  for(i in 1:n){
    v    = rep(0,2)
    v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
    v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
    v    = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
  }
  
  # Sample the weights
  w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  
  # Sample the parameters of the components
  for(k in 1:2){
    # Sample the means
    nk    = sum(cc==k)
    xsumk = sum(x[cc==k])
    tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
    mu.hat  = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
    mu[k]   = rnorm(1, mu.hat, sqrt(tau2.hat))
    # Sample the variances
    dd.star = dd + nk/2
    qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
    sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
  }
  
  # Store samples
  cc.out[s,]    = cc
  w.out[s]     = w
  mu.out[s,]    = mu
  sigma.out[s,] = sigma
  logpost[s] = 0
  for(i in 1:n){
    if(cc[i] == 1){
      logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
    }else{
      logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
    }
  }
  logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
  for(k in 1:2){
    logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
  }
  if(s/500==floor(s/500)){
    print(paste("s =",s))
  }
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx  = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
  density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) + 
                               (1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

  • 0点 No
  • 1点 Yes



2.3.2 2nd Peer

2.3.2.1 Assignment

Provide an MCMC algorithm to fit the two-component, location-and-scale mixture of Gaussians.

### Loading data and setting up global variables
library(MASS)
library(MCMCpack)
data(faithful)
KK = 2          # As asked
x  = faithful[,1]
n  = length(x)
set.seed(781209)

### Get a "Bayesian" kernel density estimator based on the same location mixture of 2 normals
## Priors set up using an "empirical Bayes" approach
aa  = rep(1,KK)  
eta = mean(x)    
tau = sqrt(var(x))
dd  = 2
qq  = var(x)/KK

w     = rep(1,KK)/KK
mu    = rnorm(KK, mean(x), sd(x))
sigma = rep(sd(x)/KK, KK)
cc    = sample(1:KK, n, replace=T, prob=w)

## Number of iterations of the sampler
rrr   = 12000
burn  = 2000

## Storing the samples
cc.out    = array(0, dim=c(rrr, n))
w.out     = array(0, dim=c(rrr, KK))
mu.out    = array(0, dim=c(rrr, KK))
sigma.out = array(0, dim=c(rrr, KK))
logpost   = rep(0, rrr)

for(s in 1:rrr){
  # Sample the indicators
  for(i in 1:n){
    v = rep(0,KK)
    for(k in 1:KK){
      v[k] = log(w[k]) + dnorm(x[i], mu[k], sigma[k], log=TRUE)  #Compute the log of the weights
    }
    v = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i] = sample(1:KK, 1, replace=TRUE, prob=v)
  }
  
  # Sample the weights
  w = as.vector(rdirichlet(1, aa + tabulate(cc, nbins=KK)))
  
  # Sample the means
  for(k in 1:KK){
    nk    = sum(cc==k)
    xsumk = sum(x[cc==k])
    tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
    mu.hat  = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
    mu[k]   = rnorm(1, mu.hat, sqrt(tau2.hat))
  }
  
  # Sample the variances
  for(k in 1:KK) {
    dd.star = dd + sum(cc==k)/2
    sumsq = 0
    for(i in 1:n) {
      if (cc[i]==k) {
        sumsq = sumsq + ((x[i]-mu[k])^2)
      }
    }
    qq.star = qq + sumsq/2
    sigma[k] = sqrt(rinvgamma(1, shape = dd.star, scale = qq.star))
  }
  
  # Store samples
  cc.out[s,]   = cc
  w.out[s,]    = w
  mu.out[s,]   = mu
  sigma.out[s,] = sigma
  for(i in 1:n){
    logpost[s] = logpost[s] + log(w[cc[i]]) + dnorm(x[i], mu[cc[i]], sigma[cc[i]], log=TRUE)
  }
  logpost[s] = logpost[s] + log(ddirichlet(w, aa))
  for(k in 1:KK){
    logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log=TRUE)
    logpost[s] = logpost[s] + log(dinvgamma(sigma[k]^2, shape = dd, scale = qq))
  }
  if(s/500==floor(s/500)){
    print(paste("s =",s, w, mu, sigma, logpost[s]))
  }
}
## [1] "s = 500 0.65972914904991 4.28469003268875 0.405452035141612 -284.691001291928" "s = 500 0.34027085095009 1.99006962679139 0.250995183146609 -284.691001291928"
## [1] "s = 1000 0.595942847037826 4.27800551958405 0.377348480705164 -286.969076980901" "s = 1000 0.404057152962174 2.03894151043637 0.30982681142742 -286.969076980901" 
## [1] "s = 1500 0.647284731246712 4.32792289508716 0.42383087208717 -283.984683957035"  "s = 1500 0.352715268753288 2.06264939688377 0.293003644635685 -283.984683957035"
## [1] "s = 2000 0.681875974292433 4.31964666156271 0.411719777867089 -284.392728508546" "s = 2000 0.318124025707567 2.05510718985465 0.295371255576304 -284.392728508546"
## [1] "s = 2500 0.601103182730055 4.24400160743064 0.446830400212041 -285.459906307434" "s = 2500 0.398896817269945 2.03587612253421 0.241473527538476 -285.459906307434"
## [1] "s = 3000 0.661235314614079 4.28417774604423 0.399650525878003 -285.992706414487" "s = 3000 0.338764685385921 2.07827156671072 0.297205073536538 -285.992706414487"
## [1] "s = 3500 0.644527213106916 4.28812055333749 0.441452138427281 -283.120726285281" "s = 3500 0.355472786893084 2.01146104149005 0.258057420679898 -283.120726285281"
## [1] "s = 4000 0.687810647492561 4.2844578377926 0.458980223432612 -287.498421765689"  "s = 4000 0.312189352507439 2.07489936081685 0.259843371397385 -287.498421765689"
## [1] "s = 4500 0.675398367857913 4.23091315341658 0.410829138479542 -285.115602530381" "s = 4500 0.324601632142087 2.00548178399042 0.263106146522126 -285.115602530381"
## [1] "s = 5000 0.623677983258687 4.2986070783669 0.460990184040201 -287.679382635602"  "s = 5000 0.376322016741313 2.04273694828362 0.296539486159871 -287.679382635602"
## [1] "s = 5500 0.650641729955276 4.30297831389609 0.410588171771523 -284.091323207027" "s = 5500 0.349358270044724 2.05603614578286 0.314091813936972 -284.091323207027"
## [1] "s = 6000 0.691002769510701 4.30695193759895 0.43643062048752 -285.25059244604"  "s = 6000 0.308997230489299 2.00724366068866 0.258042136715897 -285.25059244604"
## [1] "s = 6500 0.639653531193993 4.30228067466188 0.466022578492273 -289.071680298941" "s = 6500 0.360346468806007 2.0370347870673 0.223458414339796 -289.071680298941" 
## [1] "s = 7000 0.657266286407271 4.25015359390783 0.408697090783594 -285.609121794953" "s = 7000 0.342733713592729 1.98352245342408 0.281123375854885 -285.609121794953"
## [1] "s = 7500 0.582150245729404 4.27752884035096 0.405335835898845 -288.834567817749" "s = 7500 0.417849754270596 2.03280535197236 0.302512079952417 -288.834567817749"
## [1] "s = 8000 0.597745067402946 4.30797847743645 0.426676105386858 -286.886008869559" "s = 8000 0.402254932597054 1.97410330661178 0.266695259035665 -286.886008869559"
## [1] "s = 8500 0.67048977917052 4.34923962430531 0.440934482727733 -287.976277284725" "s = 8500 0.32951022082948 1.99450886009298 0.253385375719033 -287.976277284725"
## [1] "s = 9000 0.676526152024463 4.32975227643208 0.424164941801625 -286.420409681512" "s = 9000 0.323473847975537 2.03912276672083 0.304347017719118 -286.420409681512"
## [1] "s = 9500 0.675252958080913 4.32408885533291 0.42428978538769 -284.557683174045"  "s = 9500 0.324747041919087 2.01338154879308 0.269521567887484 -284.557683174045"
## [1] "s = 10000 0.657398814367743 4.2933157987876 0.458707092652969 -283.811172985985"  "s = 10000 0.342601185632257 2.03282618126385 0.248515730931108 -283.811172985985"
## [1] "s = 10500 0.662058917060533 4.31483111471995 0.427552521301718 -284.108178852953" "s = 10500 0.337941082939467 2.00950256527503 0.279659447389632 -284.108178852953"
## [1] "s = 11000 0.669552336247041 4.36928639963658 0.427369327130037 -288.518332324442" "s = 11000 0.330447663752959 1.97992561769347 0.296908732685882 -288.518332324442"
## [1] "s = 11500 0.627693164995977 4.20937955935309 0.437031029444333 -287.287695169602" "s = 11500 0.372306835004023 2.0025029877848 0.248798427816352 -287.287695169602" 
## [1] "s = 12000 0.706166055902923 4.32329970006521 0.423366944818915 -291.385728914099" "s = 12000 0.293833944097077 2.09369936916931 0.272464988230699 -291.385728914099"

Provide code to generate the density estimate on a grid consisting of 150 points in the interval [0,7].

Eruptions = seq(0,7, length=150)
## Compute the samples of the density over a dense grid
density.mcmc = array(0, dim=c(rrr-burn,length(Eruptions)))
for(s in 1:(rrr-burn)){
  for(k in 1:KK){
    density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn,k]*dnorm(Eruptions,mu.out[s+burn,k],sigma.out[s+burn,k])
  }
}
density.mcmc.m = apply(density.mcmc , 2, mean)

## Plot Bayesian estimate with pointwise credible bands along with kernel density estimate and frequentist point estimate
colscale = c("black", "blue", "red")
#yy = density(x)
yy = density(Eruptions)
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(Eruptions, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Duration", ylab="Density")
polygon(c(Eruptions,rev(Eruptions)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(Eruptions, density.mcmc.m, col=colscale[3], lwd=2)
#lines(xx, density.EM, col=colscale[2], lty=2, lwd=2)
#lines(yy, col=colscale[3], lty=3, lwd=2)
points(x, rep(0,n))

Provide the a graph of the density estimate on the interval [0,7].

Source :

2.3.2.2 Marking

Is the sampler for \(c\) correct?

The full conditional for \(c_i\) is given by,

\[ \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \]

(This assumes that the mixture weights are parameterized as

for(i in 1:n){
  v    = rep(0,2)
  v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
  v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
  v    = exp(v - max(v))/sum(exp(v - max(v)))
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
  • 0点 No
  • 1点 Yes

Is the sampler for \(\omega\) correct?

The full conditional is simply

\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]

The following one line of code implements it:

w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  • 0点 No
  • 1点 Yes

Are the sampler for \(\mu\) correct?

\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]

The following code implements the sampler:

# Sample the means
for(k in 1:2){
  nk       = sum(cc==k)
  xsumk    = sum(x[cc==k])
  tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
  mu.hat   = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
  mu[k]    = rnorm(1, mu.hat, sqrt(tau2.hat))
}

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Are the samplers for \(\sigma_k\) correct?

The full conditional takes the form:

\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]

The following code implements the sampler:

# Sample the variances
for(k in 1:2){
  nk    = sum(cc==k)
  dd.star = dd + nk/2
  qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
  sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Does the code to generate the density estimate appear correct?

This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

Note that this will generate not only a point estimate, but also pointwise credible intervals.

  • 0点 No
  • 1点 Yes

Does the graph appear to be correct?

The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.

The full code for the MCMC algorithm follows:

x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))

## Priors set up using an "empirical Bayes" approach
aa  = rep(1,2)  
eta = mean(x)    
tau = sqrt(var(x))
dd  = 2
qq  = var(x)/2

## Initialize the parameters
w     = 1/2
mu    = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc    = sample(1:2, n, replace=T, prob=c(1/2,1/2))

## Number of iterations of the sampler
rrr   = 12000
burn  = 2000

## Storing the samples
cc.out    = array(0, dim=c(rrr, n))
w.out     = rep(0, rrr)
mu.out    = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost   = rep(0, rrr)

for(s in 1:rrr){
  # Sample the indicators
  for(i in 1:n){
    v    = rep(0,2)
    v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
    v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
    v    = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
  }
  
  # Sample the weights
  w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  
  # Sample the parameters of the components
  for(k in 1:2){
    # Sample the means
    nk    = sum(cc==k)
    xsumk = sum(x[cc==k])
    tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
    mu.hat  = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
    mu[k]   = rnorm(1, mu.hat, sqrt(tau2.hat))
    # Sample the variances
    dd.star = dd + nk/2
    qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
    sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
  }
  
  # Store samples
  cc.out[s,]    = cc
  w.out[s]     = w
  mu.out[s,]    = mu
  sigma.out[s,] = sigma
  logpost[s] = 0
  for(i in 1:n){
    if(cc[i] == 1){
      logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
    }else{
      logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
    }
  }
  logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
  for(k in 1:2){
    logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
  }
  if(s/500==floor(s/500)){
    print(paste("s =",s))
  }
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx  = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
  density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) + 
                               (1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

  • 0点 No
  • 1点 Yes



2.3.3 3rd Peer

2.3.3.1 Assignment

2.3.3.2 Marking

Is the sampler for \(c\) correct?

The full conditional for \(c_i\) is given by,

\[ \Pr(c_i = 1|...) = \frac{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}}{\omega\frac{1}{\sqrt[]{2\pi\sigma_1}}exp\left\{-\frac{1}{2\sigma^2_1}(x_i-\mu_1)^2\right\}+(1-\omega)\frac{1}{\sqrt[]{2\pi\sigma_2}}exp\left\{-\frac{1}{2\sigma^2_2}(x_i-\mu_2)^2\right\}} \]

(This assumes that the mixture weights are parameterized as

for(i in 1:n){
  v    = rep(0,2)
  v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
  v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
  v    = exp(v - max(v))/sum(exp(v - max(v)))
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}
  • 0点 No
  • 1点 Yes

Is the sampler for \(\omega\) correct?

The full conditional is simply

\[ \omega \mid \cdots \mbox{Beta} \left(\alpha_1 + \sum_{i=1}^{n} \mathbf{1}(c_i=1), \alpha_2 + \sum_{i=1}^{n} \mathbf{1}(c_i=2)\right) \]

The following one line of code implements it:

w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  • 0点 No
  • 1点 Yes

Are the sampler for \(\mu\) correct?

\[ \mu_k \mid \cdots \mbox{Normal}\left( \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\left[ \frac{\sum_{i:c_i=k} x_i}{\sigma_k^2} + \frac{\eta}{\tau^2} \right], \left[ \frac{n_k}{\sigma_k^2} + \frac{1}{\tau^2} \right]^{-1}\right) \]

The following code implements the sampler:

# Sample the means
for(k in 1:2){
  nk       = sum(cc==k)
  xsumk    = sum(x[cc==k])
  tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
  mu.hat   = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
  mu[k]    = rnorm(1, mu.hat, sqrt(tau2.hat))
}

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Are the samplers for \(\sigma_k\) correct?

The full conditional takes the form:

\[ 1/\sigma^2_k \mid \cdots \mbox{Gamma}\left(d + \frac{n_k}{2} , q + \frac{1}{2} \sum_{i:c_k=k} (x_i - \mu_k)^2 \right) \]

The following code implements the sampler:

# Sample the variances
for(k in 1:2){
  nk    = sum(cc==k)
  dd.star = dd + nk/2
  qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
  sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
}

Remember, however, that there are various ways in which the implementation could proceed.

  • 0点 No
  • 1点 Yes

Does the code to generate the density estimate appear correct?

This is an example that, as before, assumes that the weights are parameterized as \(\omega\) and \(1- \omega\).

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

Note that this will generate not only a point estimate, but also pointwise credible intervals.

  • 0点 No
  • 1点 Yes

Does the graph appear to be correct?

The density estimate should look similar to the Figure below. Note that the variance of both components is clearly different.

The full code for the MCMC algorithm follows:

x = faithful$eruptions
n = length(x)
plot(density(x))
points(x,rep(0,n))

## Priors set up using an "empirical Bayes" approach
aa  = rep(1,2)  
eta = mean(x)    
tau = sqrt(var(x))
dd  = 2
qq  = var(x)/2

## Initialize the parameters
w     = 1/2
mu    = rnorm(2, mean(x), sd(x))
sigma = rep(sd(x)/2,2)
cc    = sample(1:2, n, replace=T, prob=c(1/2,1/2))

## Number of iterations of the sampler
rrr   = 12000
burn  = 2000

## Storing the samples
cc.out    = array(0, dim=c(rrr, n))
w.out     = rep(0, rrr)
mu.out    = array(0, dim=c(rrr, 2))
sigma.out = array(0, dim=c(rrr, 2))
logpost   = rep(0, rrr)

for(s in 1:rrr){
  # Sample the indicators
  for(i in 1:n){
    v    = rep(0,2)
    v[1] = log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)  #Compute the log of the weights
    v[2] = log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)  #Compute the log of the weights
    v    = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
  }
  
  # Sample the weights
  w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  
  # Sample the parameters of the components
  for(k in 1:2){
    # Sample the means
    nk    = sum(cc==k)
    xsumk = sum(x[cc==k])
    tau2.hat = 1/(nk/sigma[k]^2 + 1/tau^2)
    mu.hat  = tau2.hat*(xsumk/sigma[k]^2 + eta/tau^2)
    mu[k]   = rnorm(1, mu.hat, sqrt(tau2.hat))
    # Sample the variances
    dd.star = dd + nk/2
    qq.star = qq + sum((x[cc==k] - mu[k])^2)/2
    sigma[k] = sqrt(1/rgamma(1, dd.star, qq.star))
  }
  
  # Store samples
  cc.out[s,]    = cc
  w.out[s]     = w
  mu.out[s,]    = mu
  sigma.out[s,] = sigma
  logpost[s] = 0
  for(i in 1:n){
    if(cc[i] == 1){
      logpost[s] = logpost[s] + log(w) + dnorm(x[i], mu[1], sigma[1], log=TRUE)
    }else{
      logpost[s] = logpost[s] + log(1-w) + dnorm(x[i], mu[2], sigma[2], log=TRUE)
    }
  }
  logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2], log=TRUE)
  for(k in 1:2){
    logpost[s] = logpost[s] + dnorm(mu[k], eta, tau, log = T) + dgamma(1/sigma[k]^2, dd, qq)/sigma[k]^4
  }
  if(s/500==floor(s/500)){
    print(paste("s =",s))
  }
}
## [1] "s = 500"
## [1] "s = 1000"
## [1] "s = 1500"
## [1] "s = 2000"
## [1] "s = 2500"
## [1] "s = 3000"
## [1] "s = 3500"
## [1] "s = 4000"
## [1] "s = 4500"
## [1] "s = 5000"
## [1] "s = 5500"
## [1] "s = 6000"
## [1] "s = 6500"
## [1] "s = 7000"
## [1] "s = 7500"
## [1] "s = 8000"
## [1] "s = 8500"
## [1] "s = 9000"
## [1] "s = 9500"
## [1] "s = 10000"
## [1] "s = 10500"
## [1] "s = 11000"
## [1] "s = 11500"
## [1] "s = 12000"
## Compute the samples of the density over a dense grid
xx  = seq(0,7,length=150)
density.mcmc = array(0, dim=c(rrr-burn,length(xx)))
for(s in 1:(rrr-burn)){
  density.mcmc[s,] = density.mcmc[s,] + w.out[s+burn]*dnorm(xx,mu.out[s+burn,1],sigma.out[s+burn,1]) + 
                               (1-w.out[s+burn])*dnorm(xx,mu.out[s+burn,2],sigma.out[s+burn,2])
}
density.mcmc.m = apply(density.mcmc , 2, mean)

## Plot Bayesian estimate with pointwise credible bands
density.mcmc.lq = apply(density.mcmc, 2, quantile, 0.025)
density.mcmc.uq = apply(density.mcmc, 2, quantile, 0.975)
plot(xx, density.mcmc.m, type="n",ylim=c(0,max(density.mcmc.uq)),xlab="Eruptions", ylab="Density")
polygon(c(xx,rev(xx)), c(density.mcmc.lq, rev(density.mcmc.uq)), col="grey", border="grey")
lines(xx, density.mcmc.m, col="black", lwd=2)
points(x, rep(0,n))

  • 0点 No
  • 1点 Yes



2.4 ディスカッション



3 Appendix

3.1 Blooper

3.2 Documenting File Creation

It’s useful to record some information about how your file was created.

  • File creation date: 2021-05-21
  • File latest updated date: 2021-05-22
  • R version 4.1.0 (2021-05-18)
  • rmarkdown package version: 2.8
  • File version: 1.0.0
  • Author Profile: ®γσ, Eng Lian Hu
  • GitHub: Source Code
  • Additional session information:
suppressMessages(require('dplyr', quietly = TRUE))
suppressMessages(require('magrittr', quietly = TRUE))
suppressMessages(require('formattable', quietly = TRUE))
suppressMessages(require('knitr', quietly = TRUE))
suppressMessages(require('kableExtra', quietly = TRUE))

sys1 <- devtools::session_info()$platform %>% 
  unlist %>% data.frame(Category = names(.), session_info = .)
rownames(sys1) <- NULL

sys2 <- data.frame(Sys.info()) %>% 
  dplyr::mutate(Category = rownames(.)) %>% .[2:1]
names(sys2)[2] <- c('Sys.info')
rownames(sys2) <- NULL

if (nrow(sys1) == 9 & nrow(sys2) == 8) {
  sys2 %<>% rbind(., data.frame(
  Category = 'Current time', 
  Sys.info = paste(as.character(lubridate::now('Asia/Tokyo')), 'JST🗾')))
} else {
  sys1 %<>% rbind(., data.frame(
  Category = 'Current time', 
  session_info = paste(as.character(lubridate::now('Asia/Tokyo')), 'JST🗾')))
}

sys <- cbind(sys1, sys2) %>% 
  kbl(caption = 'Additional session information:') %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive')) %>% 
  row_spec(0, background = 'DimGrey', color = 'yellow') %>% 
  column_spec(1, background = 'CornflowerBlue', color = 'red') %>% 
  column_spec(2, background = 'grey', color = 'black') %>% 
  column_spec(3, background = 'CornflowerBlue', color = 'blue') %>% 
  column_spec(4, background = 'grey', color = 'white') %>% 
  row_spec(9, bold = T, color = 'yellow', background = '#D7261E')

rm(sys1, sys2)
sys
Additional session information:
Category session_info Category Sys.info
version R version 4.1.0 (2021-05-18) sysname Linux
os Ubuntu 20.04.2 LTS release 5.8.0-54-generic
system x86_64, linux-gnu version #61~20.04.1-Ubuntu SMP Thu May 13 00:05:49 UTC 2021
ui X11 nodename Scibrokes-Trading
language en machine x86_64
collate en_US.UTF-8 login englianhu
ctype en_US.UTF-8 user englianhu
tz Asia/Tokyo effective_user englianhu
date 2021-05-22 Current time 2021-05-22 04:27:25 JST🗾